home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / GAMMA.C < prev    next >
C/C++ Source or Header  |  1995-10-03  |  14KB  |  626 lines
  1. /*                            gamma.c
  2.  *
  3.  *    Gamma function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, gamma();
  10.  * extern int sgngam;
  11.  *
  12.  * y = gamma( x );
  13.  *
  14.  *
  15.  *
  16.  * DESCRIPTION:
  17.  *
  18.  * Returns gamma function of the argument.  The result is
  19.  * correctly signed, and the sign (+1 or -1) is also
  20.  * returned in a global (extern) variable named sgngam.
  21.  * This variable is also filled in by the logarithmic gamma
  22.  * function lgam().
  23.  *
  24.  * Arguments |x| <= 34 are reduced by recurrence and the function
  25.  * approximated by a rational function of degree 6/7 in the
  26.  * interval (2,3).  Large arguments are handled by Stirling's
  27.  * formula. Large negative arguments are made positive using
  28.  * a reflection formula.
  29.  *
  30.  *
  31.  * ACCURACY:
  32.  *
  33.  *                      Relative error:
  34.  * arithmetic   domain     # trials      peak         rms
  35.  *    DEC      -34, 34      10000       1.3e-16     2.5e-17
  36.  *    IEEE    -170,-33      20000       2.3e-15     3.3e-16
  37.  *    IEEE     -33,  33     20000       9.4e-16     2.2e-16
  38.  *    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
  39.  *
  40.  * Error for arguments outside the test range will be larger
  41.  * owing to error amplification by the exponential function.
  42.  *
  43.  */
  44. /* -------------------------------------------------------------- */
  45. /*    lgam()
  46.  *
  47.  *    Natural logarithm of gamma function
  48.  *
  49.  *
  50.  *
  51.  * SYNOPSIS:
  52.  *
  53.  * double x, y, lgam();
  54.  * extern int sgngam;
  55.  *
  56.  * y = lgam( x );
  57.  *
  58.  *
  59.  *
  60.  * DESCRIPTION:
  61.  *
  62.  * Returns the base e (2.718...) logarithm of the absolute
  63.  * value of the gamma function of the argument.
  64.  * The sign (+1 or -1) of the gamma function is returned in a
  65.  * global (extern) variable named sgngam.
  66.  *
  67.  * For arguments greater than 13, the logarithm of the gamma
  68.  * function is approximated by the logarithmic version of
  69.  * Stirling's formula using a polynomial approximation of
  70.  * degree 4. Arguments between -33 and +33 are reduced by
  71.  * recurrence to the interval [2,3] of a rational approximation.
  72.  * The cosecant reflection formula is employed for arguments
  73.  * less than -33.
  74.  *
  75.  * Arguments greater than MAXLGM return MAXNUM and an error
  76.  * message.  MAXLGM = 2.035093e36 for DEC
  77.  * arithmetic or 2.556348e305 for IEEE arithmetic.
  78.  *
  79.  *
  80.  *
  81.  * ACCURACY:
  82.  *
  83.  *
  84.  * arithmetic      domain        # trials     peak         rms
  85.  *    DEC     0, 3                  7000     5.2e-17     1.3e-17
  86.  *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
  87.  *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
  88.  *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
  89.  * The error criterion was relative when the function magnitude
  90.  * was greater than one but absolute when it was less than one.
  91.  *
  92.  * The following test used the relative error criterion, though
  93.  * at certain points the relative error could be much higher than
  94.  * indicated.
  95.  *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
  96.  *
  97.  */
  98. /* -------------------------------------------------------------- */
  99. /*     gamma.c        */
  100. /*    gamma function    */
  101.  
  102. /*
  103. Cephes Math Library Release 2.2:  July, 1992
  104. Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
  105. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  106. */
  107.  
  108.  
  109. #include "mconf.h"
  110.  
  111. #ifdef UNK
  112. static double P[] = {
  113.   1.60119522476751861407E-4,
  114.   1.19135147006586384913E-3,
  115.   1.04213797561761569935E-2,
  116.   4.76367800457137231464E-2,
  117.   2.07448227648435975150E-1,
  118.   4.94214826801497100753E-1,
  119.   9.99999999999999996796E-1
  120. };
  121. static double Q[] = {
  122. -2.31581873324120129819E-5,
  123.  5.39605580493303397842E-4,
  124. -4.45641913851797240494E-3,
  125.  1.18139785222060435552E-2,
  126.  3.58236398605498653373E-2,
  127. -2.34591795718243348568E-1,
  128.  7.14304917030273074085E-2,
  129.  1.00000000000000000320E0
  130. };
  131. #define MAXGAM 171.624376956302725
  132. static double LOGPI = 1.14472988584940017414;
  133. #endif
  134.  
  135. #ifdef DEC
  136. static short P[] = {
  137. 0035047,0162701,0146301,0005234,
  138. 0035634,0023437,0032065,0176530,
  139. 0036452,0137157,0047330,0122574,
  140. 0037103,0017310,0143041,0017232,
  141. 0037524,0066516,0162563,0164605,
  142. 0037775,0004671,0146237,0014222,
  143. 0040200,0000000,0000000,0000000
  144. };
  145. static short Q[] = {
  146. 0134302,0041724,0020006,0116565,
  147. 0035415,0072121,0044251,0025634,
  148. 0136222,0003447,0035205,0121114,
  149. 0036501,0107552,0154335,0104271,
  150. 0037022,0135717,0014776,0171471,
  151. 0137560,0034324,0165024,0037021,
  152. 0037222,0045046,0047151,0161213,
  153. 0040200,0000000,0000000,0000000
  154. };
  155. #define MAXGAM 34.84425627277176174
  156. static short LPI[4] = {
  157. 0040222,0103202,0043475,0006750,
  158. };
  159. #define LOGPI *(double *)LPI
  160. #endif
  161.  
  162. #ifdef IBMPC
  163. static short P[] = {
  164. 0x2153,0x3998,0xfcb8,0x3f24,
  165. 0xbfab,0xe686,0x84e3,0x3f53,
  166. 0x14b0,0xe9db,0x57cd,0x3f85,
  167. 0x23d3,0x18c4,0x63d9,0x3fa8,
  168. 0x7d31,0xdcae,0x8da9,0x3fca,
  169. 0xe312,0x3993,0xa137,0x3fdf,
  170. 0x0000,0x0000,0x0000,0x3ff0
  171. };
  172. static short Q[] = {
  173. 0xd3af,0x8400,0x487a,0xbef8,
  174. 0x2573,0x2915,0xae8a,0x3f41,
  175. 0xb44a,0xe750,0x40e4,0xbf72,
  176. 0xb117,0x5b1b,0x31ed,0x3f88,
  177. 0xde67,0xe33f,0x5779,0x3fa2,
  178. 0x87c2,0x9d42,0x071a,0xbfce,
  179. 0x3c51,0xc9cd,0x4944,0x3fb2,
  180. 0x0000,0x0000,0x0000,0x3ff0
  181. };
  182. #define MAXGAM 171.624376956302725
  183. static short LPI[4] = {
  184. 0xa1bd,0x48e7,0x50d0,0x3ff2,
  185. };
  186. #define LOGPI *(double *)LPI
  187. #endif
  188.  
  189. #ifdef MIEEE
  190. static short P[] = {
  191. 0x3f24,0xfcb8,0x3998,0x2153,
  192. 0x3f53,0x84e3,0xe686,0xbfab,
  193. 0x3f85,0x57cd,0xe9db,0x14b0,
  194. 0x3fa8,0x63d9,0x18c4,0x23d3,
  195. 0x3fca,0x8da9,0xdcae,0x7d31,
  196. 0x3fdf,0xa137,0x3993,0xe312,
  197. 0x3ff0,0x0000,0x0000,0x0000
  198. };
  199. static short Q[] = {
  200. 0xbef8,0x487a,0x8400,0xd3af,
  201. 0x3f41,0xae8a,0x2915,0x2573,
  202. 0xbf72,0x40e4,0xe750,0xb44a,
  203. 0x3f88,0x31ed,0x5b1b,0xb117,
  204. 0x3fa2,0x5779,0xe33f,0xde67,
  205. 0xbfce,0x071a,0x9d42,0x87c2,
  206. 0x3fb2,0x4944,0xc9cd,0x3c51,
  207. 0x3ff0,0x0000,0x0000,0x0000
  208. };
  209. #define MAXGAM 171.624376956302725
  210. static short LPI[4] = {
  211. 0x3ff2,0x50d0,0x48e7,0xa1bd,
  212. };
  213. #define LOGPI *(double *)LPI
  214. #endif
  215.  
  216. /* Stirling's formula for the gamma function */
  217. #if UNK
  218. static double STIR[5] = {
  219.  7.87311395793093628397E-4,
  220. -2.29549961613378126380E-4,
  221. -2.68132617805781232825E-3,
  222.  3.47222221605458667310E-3,
  223.  8.33333333333482257126E-2,
  224. };
  225. #define MAXSTIR 143.01608
  226. static double SQTPI = 2.50662827463100050242E0;
  227. #endif
  228. #if DEC
  229. static short STIR[20] = {
  230. 0035516,0061622,0144553,0112224,
  231. 0135160,0131531,0037460,0165740,
  232. 0136057,0134460,0037242,0077270,
  233. 0036143,0107070,0156306,0027751,
  234. 0037252,0125252,0125252,0146064,
  235. };
  236. #define MAXSTIR 26.77
  237. static short SQT[4] = {
  238. 0040440,0066230,0177661,0034055,
  239. };
  240. #define SQTPI *(double *)SQT
  241. #endif
  242. #if IBMPC
  243. static short STIR[20] = {
  244. 0x7293,0x592d,0xcc72,0x3f49,
  245. 0x1d7c,0x27e6,0x166b,0xbf2e,
  246. 0x4fd7,0x07d4,0xf726,0xbf65,
  247. 0xc5fd,0x1b98,0x71c7,0x3f6c,
  248. 0x5986,0x5555,0x5555,0x3fb5,
  249. };
  250. #define MAXSTIR 143.01608
  251. static short SQT[4] = {
  252. 0x2706,0x1ff6,0x0d93,0x4004,
  253. };
  254. #define SQTPI *(double *)SQT
  255. #endif
  256. #if MIEEE
  257. static short STIR[20] = {
  258. 0x3f49,0xcc72,0x592d,0x7293,
  259. 0xbf2e,0x166b,0x27e6,0x1d7c,
  260. 0xbf65,0xf726,0x07d4,0x4fd7,
  261. 0x3f6c,0x71c7,0x1b98,0xc5fd,
  262. 0x3fb5,0x5555,0x5555,0x5986,
  263. };
  264. #define MAXSTIR 143.01608
  265. static short SQT[4] = {
  266. 0x4004,0x0d93,0x1ff6,0x2706,
  267. };
  268. #define SQTPI *(double *)SQT
  269. #endif
  270.  
  271. int sgngam = 0;
  272. extern int sgngam;
  273. extern double MAXLOG, MAXNUM, PI;
  274.  
  275.  
  276. /* Gamma function computed by Stirling's formula.
  277.  * The polynomial STIR is valid for 33 <= x <= 172.
  278.  */
  279. # if defined(__STDC__) || defined(__PROTO__)
  280. static double
  281. stirf(double x)
  282. # else
  283. static double
  284. stirf(x)
  285. double    x;
  286. # endif
  287. {
  288. double y, w, v;
  289. double pow(), exp(), polevl();
  290. double log();
  291.  
  292. w = 1.0/x;
  293. w = 1.0 + w * polevl( w, STIR, 4 );
  294. y = exp(x);
  295. if( x > MAXSTIR )
  296.     { /* Avoid overflow in pow() */
  297.     v = pow( x, 0.5 * x - 0.25 );
  298.     y = v * (v / y);
  299.     }
  300. else
  301.     {
  302.     y = pow( x, x - 0.5 ) / y;
  303.     }
  304. y = SQTPI * y * w;
  305. return( y );
  306. }
  307.  
  308.  
  309. # if defined(__STDC__) || defined(__PROTO__)
  310. double
  311. gamma(double x)
  312. # else
  313. double
  314. gamma(x)
  315. double    x;
  316. # endif
  317. {
  318. double p, q, z;
  319. int i;
  320. double fabs(), lgam(), log(), exp(), gamma(), sin();
  321. double floor(), polevl(), p1evl(), stirf();
  322.  
  323. sgngam = 1;
  324. q = fabs(x);
  325.  
  326. if( q > 33.0 )
  327.     {
  328.     if( x < 0.0 )
  329.         {
  330.         p = floor(q);
  331.         if( p == q )
  332.             goto goverf;
  333.         i = (int)p;
  334.         if( (i & 1) == 0 )
  335.             sgngam = -1;
  336.         z = q - p;
  337.         if( z > 0.5 )
  338.             {
  339.             p += 1.0;
  340.             z = q - p;
  341.             }
  342.         z = q * sin( PI * z );
  343.         if( z == 0.0 )
  344.             {
  345. goverf:
  346.             mtherr( "gamma", OVERFLOW );
  347.             return( sgngam * MAXNUM);
  348.             }
  349.         z = fabs(z);
  350.         z = PI/(z * stirf(q) );
  351.         }
  352.     else
  353.         {
  354.         z = stirf(x);
  355.         }
  356.     return( sgngam * z );
  357.     }
  358.  
  359. z = 1.0;
  360. while( x >= 3.0 )
  361.     {
  362.     x -= 1.0;
  363.     z *= x;
  364.     }
  365.  
  366. while( x < 0.0 )
  367.     {
  368.     if( x > -1.E-9 )
  369.         goto small;
  370.     z /=x;
  371.     x += 1.0;
  372.     }
  373.  
  374. while( x < 2.0 )
  375.     {
  376.     if( x < 1.e-9 )
  377.         goto small;
  378.     z /=x;
  379.     x += 1.0;
  380.     }
  381.  
  382. if( (x == 2.0) || (x == 3.0) )
  383.     return(z);
  384.  
  385. x -= 2.0;
  386. p = polevl( x, P, 6 );
  387. q = polevl( x, Q, 7 );
  388. return( z * p / q );
  389.  
  390. small:
  391. if( x == 0.0 )
  392.     {
  393.     mtherr( "gamma", SING );
  394.     return( MAXNUM );
  395.     }
  396. else
  397.     return( z/((1.0 + 0.5772156649015329 * x) * x) );
  398. }
  399.  
  400.  
  401.  
  402. /* A[]: Stirling's formula expansion of log gamma
  403.  * B[], C[]: log gamma function between 2 and 3
  404.  */
  405. #ifdef UNK
  406. static double A[] = {
  407.  8.11614167470508450300E-4,
  408. -5.95061904284301438324E-4,
  409.  7.93650340457716943945E-4,
  410. -2.77777777730099687205E-3,
  411.  8.33333333333331927722E-2
  412. };
  413. static double B[] = {
  414. -1.37825152569120859100E3,
  415. -3.88016315134637840924E4,
  416. -3.31612992738871184744E5,
  417. -1.16237097492762307383E6,
  418. -1.72173700820839662146E6,
  419. -8.53555664245765465627E5
  420. };
  421. static double C[] = {
  422.  1.00000000000000000000E0,
  423. -3.51815701436523470549E2,
  424. -1.70642106651881159223E4,
  425. -2.20528590553854454839E5,
  426. -1.13933444367982507207E6,
  427. -2.53252307177582951285E6,
  428. -2.01889141433532773231E6
  429. };
  430. /* log( sqrt( 2*pi ) ) */
  431. static double LS2PI  =  0.91893853320467274178;
  432. #define MAXLGM 2.556348e305
  433. #endif
  434.  
  435. #ifdef DEC
  436. static short A[] = {
  437. 0035524,0141201,0034633,0031405,
  438. 0135433,0176755,0126007,0045030,
  439. 0035520,0006371,0003342,0172730,
  440. 0136066,0005540,0132605,0026407,
  441. 0037252,0125252,0125252,0125132
  442. };
  443. static short B[] = {
  444. 0142654,0044014,0077633,0035410,
  445. 0144027,0110641,0125335,0144760,
  446. 0144641,0165637,0142204,0047447,
  447. 0145215,0162027,0146246,0155211,
  448. 0145322,0026110,0010317,0110130,
  449. 0145120,0061472,0120300,0025363
  450. };
  451. static short C[] = {
  452. /*0040200,0000000,0000000,0000000*/
  453. 0142257,0164150,0163630,0112622,
  454. 0143605,0050153,0156116,0135272,
  455. 0144527,0056045,0145642,0062332,
  456. 0145213,0012063,0106250,0001025,
  457. 0145432,0111254,0044577,0115142,
  458. 0145366,0071133,0050217,0005122
  459. };
  460. /* log( sqrt( 2*pi ) ) */
  461. static short LS2P[] = {040153,037616,041445,0172645,};
  462. #define LS2PI *(double *)LS2P
  463. #define MAXLGM 2.035093e36
  464. #endif
  465.  
  466. #ifdef IBMPC
  467. static short A[] = {
  468. 0x6661,0x2733,0x9850,0x3f4a,
  469. 0xe943,0xb580,0x7fbd,0xbf43,
  470. 0x5ebb,0x20dc,0x019f,0x3f4a,
  471. 0xa5a1,0x16b0,0xc16c,0xbf66,
  472. 0x554b,0x5555,0x5555,0x3fb5
  473. };
  474. static short B[] = {
  475. 0x6761,0x8ff3,0x8901,0xc095,
  476. 0xb93e,0x355b,0xf234,0xc0e2,
  477. 0x89e5,0xf890,0x3d73,0xc114,
  478. 0xdb51,0xf994,0xbc82,0xc131,
  479. 0xf20b,0x0219,0x4589,0xc13a,
  480. 0x055e,0x5418,0x0c67,0xc12a
  481. };
  482. static short C[] = {
  483. /*0x0000,0x0000,0x0000,0x3ff0,*/
  484. 0x12b2,0x1cf3,0xfd0d,0xc075,
  485. 0xd757,0x7b89,0xaa0d,0xc0d0,
  486. 0x4c9b,0xb974,0xeb84,0xc10a,
  487. 0x0043,0x7195,0x6286,0xc131,
  488. 0xf34c,0x892f,0x5255,0xc143,
  489. 0xe14a,0x6a11,0xce4b,0xc13e
  490. };
  491. /* log( sqrt( 2*pi ) ) */
  492. static short LS2P[] = {
  493. 0xbeb5,0xc864,0x67f1,0x3fed
  494. };
  495. #define LS2PI *(double *)LS2P
  496. #define MAXLGM 2.556348e305
  497. #endif
  498.  
  499. #ifdef MIEEE
  500. static short A[] = {
  501. 0x3f4a,0x9850,0x2733,0x6661,
  502. 0xbf43,0x7fbd,0xb580,0xe943,
  503. 0x3f4a,0x019f,0x20dc,0x5ebb,
  504. 0xbf66,0xc16c,0x16b0,0xa5a1,
  505. 0x3fb5,0x5555,0x5555,0x554b
  506. };
  507. static short B[] = {
  508. 0xc095,0x8901,0x8ff3,0x6761,
  509. 0xc0e2,0xf234,0x355b,0xb93e,
  510. 0xc114,0x3d73,0xf890,0x89e5,
  511. 0xc131,0xbc82,0xf994,0xdb51,
  512. 0xc13a,0x4589,0x0219,0xf20b,
  513. 0xc12a,0x0c67,0x5418,0x055e
  514. };
  515. static short C[] = {
  516. 0xc075,0xfd0d,0x1cf3,0x12b2,
  517. 0xc0d0,0xaa0d,0x7b89,0xd757,
  518. 0xc10a,0xeb84,0xb974,0x4c9b,
  519. 0xc131,0x6286,0x7195,0x0043,
  520. 0xc143,0x5255,0x892f,0xf34c,
  521. 0xc13e,0xce4b,0x6a11,0xe14a
  522. };
  523. /* log( sqrt( 2*pi ) ) */
  524. static short LS2P[] = {
  525. 0x3fed,0x67f1,0xc864,0xbeb5
  526. };
  527. #define LS2PI *(double *)LS2P
  528. #define MAXLGM 2.556348e305
  529. #endif
  530.  
  531.  
  532. /* Logarithm of gamma function */
  533.  
  534.  
  535. # if defined(__STDC__) || defined(__PROTO__)
  536. double
  537. lgam(double x)
  538. # else
  539. double
  540. lgam(x)
  541. double    x;
  542. # endif
  543. {
  544. double p, q, w, z;
  545. int i;
  546. double fabs(), sin(), log(), gamma(), polevl();
  547. double p1evl(), floor();
  548.  
  549. sgngam = 1;
  550.  
  551. if( x < -34.0 )
  552.     {
  553.     q = -x;
  554.     w = lgam(q); /* note this modifies sgngam! */
  555.     p = floor(q);
  556.     if( p == q )
  557.         goto loverf;
  558.     i = (int)p;
  559.     if( (i & 1) == 0 )
  560.         sgngam = -1;
  561.     else
  562.         sgngam = 1;
  563.     z = q - p;
  564.     if( z > 0.5 )
  565.         {
  566.         p += 1.0;
  567.         z = p - q;
  568.         }
  569.     z = q * sin( PI * z );
  570.     if( z == 0.0 )
  571.         goto loverf;
  572. /*    z = log(PI) - log( z ) - w;*/
  573.     z = LOGPI - log( z ) - w;
  574.     return( z );
  575.     }
  576.  
  577. if( x < 13.0 )
  578.     {
  579.     z = 1.0;
  580.     while( x >= 3.0 )
  581.         {
  582.         x -= 1.0;
  583.         z *= x;
  584.         }
  585.     while( x < 2.0 )
  586.         {
  587.         if( x == 0.0 )
  588.             goto loverf;
  589.         z /=x;
  590.         x += 1.0;
  591.         }
  592.     if( z < 0.0 )
  593.         {
  594.         sgngam = -1;
  595.         z = -z;
  596.         }
  597.     else
  598.         sgngam = 1;
  599.     if( x == 2.0 )
  600.         return( log(z) );
  601.     x -= 2.0;
  602.     p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
  603.     return( log(z) + p );
  604.     }
  605.  
  606. if( x > MAXLGM )
  607.     {
  608. loverf:
  609.     mtherr( "lgam", OVERFLOW );
  610.     return( sgngam * MAXNUM );
  611.     }
  612.  
  613. q = ( x - 0.5 ) * log(x) - x + LS2PI;
  614. if( x > 1.0e8 )
  615.     return( q );
  616.  
  617. p = 1.0/(x*x);
  618. if( x >= 1000.0 )
  619.     q += ((   7.9365079365079365079365e-4 * p
  620.         - 2.7777777777777777777778e-3) *p
  621.         + 0.0833333333333333333333) / x;
  622. else
  623.     q += polevl( p, A, 4 ) / x;
  624. return( q );
  625. }
  626.